c1.est=function(data,y,m,M=NULL,N,param){
  data=as.name(data); N=N; y=y; m=m; n=length(y)
  fpc=(1-n/N); ybar=sum(y)/sum(m)
  s2r=sum((y-ybar*m)^2)/(n-1)
  plot(m,y,main=data)
  if((!is.null(M))){
    mbar=M/N
    if(param=='mean'){
      est=ybar
      vhat=fpc*(s2r/(n*mbar^2))
      bound=2*sqrt(vhat)
      lower=est-bound; upper=est+bound
    } else if(param=='total'){
      est=M*ybar
      vhat=M^2*fpc*(s2r/(n*mbar^2))
      bound=2*sqrt(vhat)
      lower=est-bound; upper=est+bound
    } else{
      est=sum(y)/sum(m)
      s2p=sum((y-m*est)^2)/(n-1)
      vhat=fpc*(s2p/(n*mbar^2))
      bound=2*sqrt(vhat)
      l=est-bound; u=est+bound
      lower=sprintf("%.2f%%",l*100)
      upper=sprintf("%.2f%%",u*100)
    }
  } else{
    mbar=sum(m)/n
    ybart=sum(y)/n
    s2t=sum((y-ybart)^2)/(n-1)
    if(param=='mean'){
      est=ybart
      vhat=fpc*(s2t/(n*mbar^2))
      bound=2*sqrt(vhat)
      lower=est-bound; upper=est+bound
    } else if(param=='total'){
      est=N*ybart
      vhat=N^2*fpc*(s2t/n)
      bound=2*sqrt(vhat)
      lower=est-bound; upper=est+bound
    } else{
      est=sum(y)/sum(m)
      s2p=sum((y-m*est)^2)/(n-1)
      vhat=fpc*(s2p/(n*mbar^2))
      bound=2*sqrt(vhat)
      l=est-bound; u=est+bound
      lower=sprintf("%.2f%%",l*100)
      upper=sprintf("%.2f%%",u*100)
    }
  } 
  bias=fpc*((1/n)*(sd(m)/mean(m))^2-cor(m,y)*(sd(m)/mean(m))*(sd(y)/mean(y)))
  cat("","\n","Results from 1-stage Cluster sample","\n","Data =",data,"\n","\n","N =",N,"n =",n,
      "FPC =",fpc,"\n","Relative bias =",bias,"\n","Estimate of",param,"=",est,"\n",
      "Vhat(",param,") =",vhat,"\n","Bound =",bound,"\n","Lower Bound =",lower,"Upper Bound =",
      upper,"\n","")
  results=list(est=est,data=data,n=n,N=N,fpc=fpc,vhat=vhat,bound=bound,lower=lower,upper=upper,
               param=param,bias=bias)
}
# to use the function with its call:
# c1.est(data,y,m,M,N,param)

# data: dataset name, in quotes
# y: numeric vector
# m: cluster elements
# M: population cluster elements (single number or NULL)
# N: population sizes (single number)
# param: c('mean','total','proportion')